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ABSTRACT 

Ocean  (T,  S )  data  analysis/assimilation,  conducted  in  the  three-dimensional  physical  space,  is  a  generalized 
average  of  purely  observed  data  (data  analysis)  or  of  modeled/observed  data  (data  assimilation).  Because  of 
the  high  nonlinearity  of  the  equation  of  the  state  of  the  seawater  and  nonuniform  vertical  distribution  of  the 
observational  profile  data,  false  static  instability  may  be  generated.  A  new  analytical  conserved  adjustment 
scheme  has  been  developed  on  the  base  of  conservation  of  heat,  salt,  and  static  stability  for  the  whole  water 
column  with  predetermined  ( T ,  S)  adjustment  ratios.  A  set  of  well-posed  combined  linear  and  nonlinear 
algebraic  equations  has  been  established  and  is  solved  using  Newton’s  method.  This  new  scheme  can  be  used 
for  ocean  hydrographic  data  analysis  and  data  assimilation. 


1.  Introduction 

Raw  and  averaged  observational  hydrographic  data 
contain  substantial  regions  with  vertical  density  in¬ 
versions.  For  example,  Jackett  and  McDougall  (1995) 
found  that  the  annually  averaged  field  of  the  ocean  atlas 
of  Levitus  (1982)  had  more  than  44%  of  the  casts  pos¬ 
sessing  static  instability  at  least  at  one  level.  Flere,  the 
word  “cast”’  is  used  to  denote  a  pair  of  vertical  tem¬ 
perature  and  salinity  profiles.  A  widely  used  concept  for 
static  stability  E  is  defined  by  Lynn  and  Reid  (1968)  as 
“the  individual  density  gradient  by  vertical  displace¬ 
ment  of  a  water  parcel  (as  opposed  to  the  geometric 
density  gradient).”  For  discrete  samples  (If.  Sk)  at 
depth  zk,  k  =  1,2 , ...  ,K{k  increasing  downward),  the 
density  difference  between  two  adjacent  levels  is  taken 
after  one  is  adiabatically  displaced  to  the  depth  of  the 
other.  Computationally,  Ek  is  calculated  by 

=  P($k+ v  Tk+vzk)  ~  P($k’  Tk,  zk), 
k=l,2,  ...  ,K—1,  (1) 

where  p(Sk+i,  Tk+i,  zk )  is  the  local  potential  density  of 
the  lower  of  the  two  adjacent  levels  between  zk  and  zk+i, 
with  respect  to  the  upper  of  the  two  adjacent  levels  zk, 
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and  p  is  the  in  situ  density  to  the  depth  of  the  upper  of 
the  two  adjacent  levels  zk.  The  density  inversion  is  de¬ 
fined  by  the  occurrence  of  the  negative  value  of  Ek.  The 
minimum  static  stability  is  represented  by  Ek  =  Emm.  It 
is  not  always  possible  to  reach  zero  exactly  due  to  the 
precision  limitations  of  the  temperature  and  salinity 
values  used  (Locarnini  et  al.  2006).  As  a  result,  the 
minimum  value  for  the  static  stability  is  given  by 

E>E.,  k  =  1,2,  ...  ,K,  (2) 

where  Em jn  is  the  reference  value  for  the  minimum  static 
stability,  which  is  user-defined.  If  static  instability  occurs 
in  an  observed  or  averaged  hydrographic  cast  [i.e.,  (2)  is 
not  satisfied],  this  profile  needs  to  be  adjusted. 

The  National  Oceanographic  Data  Center  (NODC) 
uses  a  local  interactive  (T,  S)  separated  adjustment  method 
(Locarnini  et  al.  2006),  which  is  based  on  the  method 
proposed  by  Jackett  and  McDougall  (1995)  with  some 
modifications,  to  minimally  adjust  unstable  temperature 
and  salinity  profiles  (hereafter  referred  to  as  the  MA 
method) 

x  —  (Tv  T2,  . . .  ,  Tk,  SvS2,  . . .  ,SK ) 
into  stable  profiles 

x  +  Ax  =  (Tj  +  A Tp  T2  +  AT2,  . . .  ,  TK  +  A TK, 

S)  +  AAp  S2  +  A S2,  ...  ,  SK  +  A SK). 
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After  assuming  T  and  S  are  linear,  the  adjustment  is  to 
solve  the  problem: 

Minimize 1 1 Ax |j  subject  to  A  ■  (x  +  Ax)  >  7smin,  (3) 

where  the  finite-difference  approximation  of  stability  Ek 
becomes  the  inner  product  of  the  matrix  A  and  the 
profile  vector  x  +  Ax.  Obviously,  matrix  A  depends  on 
the  solution  Ax  to  the  minimization  problem  (3),  im¬ 
plying  that  the  constraints  in  (3)  are  nonlinear.  Usually, 
an  iteration  method  is  used. 

Before  deciding  which  level  to  change,  the  values  of 
dTIdz  and  dSIdz,  the  gradients  of  temperature  and  salinity 
between  two  adjacent  levels  involved  in  the  instability, 
are  examined.  This  helps  determine  if  the  temperature  or 
salinity  profile,  or  both,  are  to  be  changed  to  stabilize  the 
density  held.  If  dTIdz  <  0,  dS/dz  <  0,  only  temperature  is 
changed;  if  dTIdz  >  0,  dSIdz  >  0,  only  salinity  is  changed; 
and  it  dTIdz  <  0,  dS/dz  >  0,  both  temperature  and  salinity 
fields  are  adjusted  with  a  local  linear  trend  test  (Locarnini 
et  al.  2006).  Here,  the  z  axis  points  upward.  The  prin¬ 
ciple  is  to  stabilize  the  hydrographic  profiles  with  mini¬ 
mum  adjustment. 

The  benefit  of  using  the  MA  method  can  be  easily 
identified  from  comparison  between  two  ocean  atlases: 
the  ocean  atlas  of  Levitus  (1982;  without  MA)  and  the 
World  Ocean  Adas  2005  (Locarnini  et  al.  2006;  with 
MA).  Both  atlases  consist  of  annually  and  monthly  av¬ 
eraged  vertical  profiles  of  temperature  and  salinity  on 
a  global  1°  X  1°  grid  at  33  vertical  levels.  The  ocean  atlas 
of  Levitus  (1982)  has  considerable  casts  possessing  static 
instability;  however,  the  World  Ocean  Adas  2005  con¬ 
tains  no  profile  possessing  static  instability. 

To  eliminate  the  static  instability,  the  MA  method  does 
not  require  the  conservation  of  heat  and  salt.  Because  one 
of  the  ocean’s  important  roles  in  the  earth’s  climate  is 
heat  transport,  an  adjustment  made  without  taking  heat 
conservation  into  account  may  lead  to  errors  in  estimat¬ 
ing  the  ocean’s  impact  on  global  climate  change.  In  this 
study,  a  new  conserved  scheme  is  developed  to  simulta¬ 
neously  adjust  the  temperature  and  salinity  profiles  from 
(7*,  Sk)  to  (Tk  +  ATk,  Sk  +  AS*).  A  set  of  2 K  algebraic 
(linear  and  nonlinear)  equations  are  established  to  get 
(A  7).,  A SK)  on  the  base  of  heat  and  salt  conservation, 
predetermined  (ATk/ASK)  ratios  (or  called  adjustment 
ratios)  for  all  levels,  and  the  removal  of  static  instability 
by  adjusting  Ek  to  Ek  +  AEk  with  a  combined  conserva¬ 
tion  and  nonuniform  increment  treatment. 


this  example  is  the  1°  latitude-longitude  box  centered 
at  53.5°S,  171. 5°E  from  Levitus  et  al.  (1998).  This  is  on 
the  New  Zealand  Plateau,  with  a  bottom  depth  below 
1000  m  and  above  1100  m.  The  month  is  October,  dur¬ 
ing  the  early  austral  summer.  There  is  no  temperature  or 
salinity  data  within  the  chosen  1°  box.  Thus,  the  objec¬ 
tively  analyzed  values  in  this  1°  box  will  be  dependent  on 
the  seasonal  objectively  analyzed  held  and  the  data  in 
nearby  1°  grid  boxes.  There  is  much  more  temperature 
data  than  salinity  data  on  the  New  Zealand  Plateau  for 
October.  This  contributes  to  six  small  (on  the  order  of 
10~2  kg  m~3)  inversions  in  the  local  potential  density 
held  calculated  from  objectively  analyzed  temperature 
and  salinity  helds  (Table  1).  After  using  the  MA  method, 
the  original  and  adjusted  prohles  {Tk,  Sk,  k  =  1,  2,  ...  , 
K}  are  as  shown  in  Fig.  1,  and  the  adjusted  temperature 
and  salinity  prohles  are  listed  in  Table  2.  Readers  are 
referred  to  appendix  B  of  Locarnini  et  al.  (2006)  for 
detailed  information  on  the  stabilization  procedures. 
The  relative  root -mean  adjustment  (RRMA)  using  the 
MA  method  can  be  represented  by 


RRMA  = 


sS<4JA 


X(AS,)2 


max(7’jt)  -  min(7’/f) 
=  0.0712. 


ma  x(Sk)  —  min(AA.) 

(4) 


RRMA  represents  the  mean  adjustment  relative  to  the 
range  of  a  prohle.  The  total  heat  and  salt  changes  of  the 
water  column  within  this  1°  X  1°  grid  box  are  estimated  by 

fo  fo 


A  Q  =  Ap0cp 


AT  dz,  A(salt)=A  ASdz, 

-H  J-H 


where  p0  (=1028  kg  m~3)  is  the  characteristic  density, 
cp  (=4002  J  kg-1  K_1)  is  the  specihc  heat  for  the  sea¬ 
water,  7/  =  1000  m,  and  A  is  the  area  of  the  grid  box, 

A=(iii*)2cos<p’ 


where  R  (=6370  km)  is  the  earth’s  radius,  and  cp  (=53.5°) 
is  the  latitude  of  the  grid  box.  The  temperature  and  sa¬ 
linity  adjustments  (A 7)  A, S')  are  obtained  by  comparison 
between  Tables  1  and  2,  the  heat  and  salt  changes  of  the 
water  column  for  this  grid  box  are  calculated  by 

A Q  =  -7.0411  X  1017  J,  A(salt)  =  -0.5443  X  1010kg. 


2.  Unconserved  adjustment 

Because  one  of  the  ocean’s  important  roles  in  the  earth’s 
An  example  as  described  in  appendix  B  of  Locarnini  climate  is  transporting  heat  from  low  to  high  latitudes, 
et  al.  (2006)  is  used  for  illustration.  The  area  chosen  for  nontrivial  heat  and  salt  losses  show  that  the  unconserved 
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Table  1.  Grid  box  171.5°E,  53.5°S  Levitus  et  al.  (1998)  profiles  before  stabilization  (from  Locarnini  et  al.  2006,  Table  Bl).  Here,  the 

asterisks  in  the  last  column  indicate  the  static  instability. 


k 

Depth  (m) 

T  (°C) 

S  (ppt) 

p(Sk+ i,  Tk+ 1,  zk)  (kg  m  3) 

p(Sk,  Tk,  zf)  (kg  m  3) 

Ek  (kg  m  3) 

1 

0 

7.1667 

34.4243 

26.9476 

26.9423 

0.0054 

2 

10 

7.1489 

34.4278 

26.8982 

26.9939 

-0.0957* 

3 

20 

7.0465 

34.2880 

26.9529 

26.9443 

0.0085 

4 

30 

7.0050 

34.2914 

27.0104 

26.9990 

0.0114 

5 

50 

6.9686 

34.2991 

27.0967 

27.1028 

-0.0061* 

6 

75 

7.0604 

34.3073 

27.2406 

27.2120 

0.0286 

7 

100 

6.9753 

34.3280 

27.3892 

27.3560 

0.0332 

8 

125 

6.9218 

34.3604 

27.5164 

27.5046 

0.0117 

9 

150 

6.8919 

34.3697 

27.6000 

27.6316 

-0.0316* 

10 

200 

6.9363 

34.3364 

27.8123 

27.8302 

-0.0179* 

11 

250 

7.0962 

34.3415 

28.0295 

28.0421 

-0.0126* 

12 

300 

7.1622 

34.3367 

28.2684 

28.2593 

0.0092 

13 

400 

6.8275 

34.2852 

28.6664 

28.7281 

-0.0618* 

14 

500 

7.4001 

34.3123 

29.3699 

29.1238 

0.2461 

15 

600 

6.2133 

34.4022 

29.9386 

29.8292 

0.1094 

16 

700 

5.9186 

34.4868 

30.5869 

30.3978 

0.1891 

17 

800 

4.5426 

34.4904 

31.0754 

31.0488 

0.0266 

18 

900 

4.1263 

34.4558 

31.6539 

31.5377 

0.1162 

19 

1000 

3.3112 

34.4755 

32.1176 

adjustment  may  change  heat  transport  and  in  turn  affect 
the  overturning  thermohaline  circulation. 


£*  =  £., 
k.  min’ 


(5) 


3.  Stabilization 

The  stabilization  process  is  divided  into  three  parts:  1) 
stability  increasing  at  unstable  levels,  2)  stability  de¬ 
creasing  at  stable  levels,  and  3)  normalization  for  con¬ 
servation  of  stability  for  the  cast.  Let  static  instability 
occur  at  level  ku  k2,  . . .  ,  [i.e.,  satisfies  the  inequality 

(2)],  the  static  stability  Ek  is  increased  to  its  marginal 
stability  value  (E£), 


that  is,  the  minimal  adjustment  with  increment  of 


^k=Emm~Ek. 


Such  an  increase  of  stability  will  be  compensated  by  the 
decrease  of  stability  at  neighboring  levels  k,  ±  m  (m  =  1, 
2,  . . .)  with  skipping  the  unstable  levels  until  reaching 
the  top  and  bottom  of  the  profile. 


f£,+  -AE./2m+1 

r-.;.  I  k.±m  k- 

Eki±m  =  \E 

I  min 


if  E.+  —  A.E.  /2m+1  >  E  . 

k-±m  k-  min 


if  E,. 


—  AE.  I2m+1  <  E  .  ' 

k-  min 


The  static  stabilities  for  the  whole  profile  before  and 
after  the  adjustment  are  calculated  by 

K  K 

/=  IX,  /*=  Si*  (6) 

k=l  k= 1 


The  normalization  process  is  conducted  by 


J7**  — 
^k 


(7) 


to  keep  the  conservation  of  the  static  stability  for  each 
profile.  After  three  stabilization  processes,  the  static 
stability  is  represented  by  [see  Eq.  (1)] 


P(^k+ 1  +  ^fc+l>  Ek+\  +  ^k+V  Zk )  P($k  +  ^k’ 

Tk  +  ATk,zk)  =  E**,  k=h2,  ...  ,K-1.  (8) 

When  Emin  is  specified  [see  Eq.  (5)],  the  right-hand  side 
of  (8)  (i.e.,  /:**)  is  the  known  adjustment,  which  is  cal¬ 
culated  through  (5)— (7).  Equation  (8)  is  used  to  de¬ 
termine  the  temperature  and  salinity  adjustments  at  each 
depth  for  given  E'£*.  The  difference  between  (7)  (i.e., 
the  direct  determination  of  E'f*)  and  (8)  is  that  (7) 
shows  the  minimal  density  adjustment  to  remove  static 
instability,  and  (8)  is  to  calculate  the  ( T ,  S )  adjustment 
at  each  depth. 
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Temperature  Salinity  Identifying  instabilities 

Fig.  1.  Original  (dashed)  and  adjusted  (solid)  profiles  temperature  of  Tk,  salinity  Skl  and  static  stability  Ek  at  the  grid  box  53.5°S,  171.5°E 

using  the  MA  method  (Locarnini  et  al.  2006). 


4.  Constraints  for  temperature  and  salinity 
adjustment 

Conservation  of  heat  and  salt  for  the  adjustment  can 
be  represented  by 


•o 

ATdz  =  0, 

-h 


•0 

ASdz  =  Q, 


-h 


which  can  be  discretized  by 


(9) 


Table  2.  Grid  box  53.5°S,  171.5°E  improved  Levitus  et  al.  (1998)  profiles  after  stabilization  using  the  MA  method  (from  Locarnini  et  al. 

2006,  Table  B2). 


k 

Depth  (m) 

T  (°C) 

S  (ppt) 

p(Sk+ 1,  Tk+1,  zk)  (kg  m  3) 

p(Sk,  Tk,  zk)  (kg  m  3) 

Ek  (kg  m"3) 

1 

0 

7.1667 

34.3096 

26.8521 

26.8519 

0.0002 

2 

10 

7.1489 

34.3063 

26.8982 

26.8982 

0.0000 

3 

20 

7.0465 

34.2880 

26.9529 

26.9443 

0.0085 

4 

30 

7.0050 

34.2914 

27.0042 

26.9990 

0.0051 

5 

50 

7.0132 

34.2991 

27.0967 

27.0967 

0.0000 

6 

75 

7.0604 

34.3073 

27.2361 

27.2120 

0.0240 

7 

100 

6.9796 

34.3228 

27.3513 

27.3513 

0.0000 

8 

125 

6.9897 

34.3243 

27.4667 

27.4667 

0.0000 

9 

150 

7.0242 

34.3301 

27.5820 

27.5820 

0.0000 

10 

200 

7.0628 

34.3364 

27.8123 

27.8123 

0.0000 

11 

250 

7.0962 

34.3415 

28.0422 

28.0421 

0.0000 

12 

300 

7.0748 

34.3367 

28.2719 

28.2719 

0.0001 

13 

400 

6.8275 

34.2894 

28.7314 

28.7314 

0.0000 

14 

500 

6.9604 

34.3123 

29.3699 

29.1899 

0.1799 

15 

600 

6.2133 

34.4022 

29.9386 

29.8292 

0.1094 

16 

700 

5.9186 

34.4868 

30.5869 

30.3978 

0.1891 

17 

800 

4.5426 

34.4904 

31.0754 

31.0488 

0.0266 

18 

900 

4.1263 

34.4558 

31.6539 

31.5377 

0.1162 

19 

1000 

3.3112 

34.4755 

32.1176 
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TABLE  3.  Change  of  (Tk 

+  ATk)  (°C)  at  each  iteration  using  the  Newton's  method.  It  is  noted  that  the  iteration  converges  at  the 

third  iteration. 

k 

Depth  (m) 

/  =  o 

7  =  1 

7  =  2 

7  =  3 

7  =  4 

1 

0 

7.166  700 

7.212  634 

7.212  833 

7.212  833 

7.212  833 

2 

10 

7.148  900 

7.289  401 

7.289  072 

7.289  072 

7.289  072 

3 

20 

7.046  500 

6.818  173 

6.816  828 

6.816  828 

6.816  828 

4 

30 

7.005  000 

6.872  865 

6.872  591 

6.872  591 

6.872  591 

5 

50 

6.968  600 

6.888  794 

6.888  861 

6.888  861 

6.888  861 

6 

75 

7.060  400 

7.023  494 

7.023  712 

7.023  712 

7.023  712 

7 

100 

6.975  300 

6.977  379 

6.977  638 

6.977  638 

6.977  638 

8 

125 

6.921  800 

6.965  175 

6.965  378 

6.965  378 

6.965  378 

9 

150 

6.891  900 

6.983  992 

6.983  997 

6.983  997 

6.983  997 

10 

200 

6.936  300 

6.959  537 

6.959  779 

6.959  779 

6.959  779 

11 

250 

7.096  200 

7.125  999 

7.126  229 

7.126  229 

7.126  229 

12 

300 

7.162  200 

7.228  075 

7.228  205 

7.228  205 

7.228  205 

13 

400 

6.827  500 

6.995  044 

6.994  489 

6.994  488 

6.994  488 

14 

500 

7.400  100 

7.229  221 

7.228  652 

7.228  652 

7.228  652 

15 

600 

6.213  300 

6.129  374 

6.129  400 

6.129  400 

6.129  400 

16 

700 

5.918  600 

5.883  923 

5.884  121 

5.884  121 

5.884  121 

17 

800 

4.542  600 

4.542  873 

4.543  127 

4.543  127 

4.543  127 

18 

900 

4.126  300 

4.153  784 

4.154  020 

4.154  020 

4.154  020 

19 

1000 

3.311  200 

3.362  894 

3.363075 

3.363  075 

3.363  075 

K-l, 


K-l  , 


Y(^  +  ^+1), 


k= 1 


(Z*-W  =  0- 


(10) 

(ii) 


Equations  (10)  and  (11)  can  be  rewritten  by 


K 

'5!akATk=0' 

k= 1 


K 

^akASk  =  0, 

k=l 


where 


(12) 


1 


Zi  Z2  _  Zl  Z3 

2(Z!-Zj’  fl2“2(z1-z/f)’  ^ 


Z2  Z4 

2(zi-^)’"'’ 


K-l 


ZK-2  ZK 

2(zi  “  Zk)  ' 


ZA-1  ZA 

2(zi“za)' 


(13) 


Obviously,  we  have 

A 

/ .a.  =  1 .  a,  >0  for  k  =  l,2,...,K.  (14) 

k=  1 

The  adjustment  ratios  are  used  for  TV  —  1  levels, 

+  ykASk  =  0,  *  =  1,2,  ...  ,K-1.  (15) 

Because  temperature  and  salinity  corrections  affect  the 
density  differently,  that  is,  the  increase  (decrease)  of 
temperature  (salinity)  decreases  (increases)  the  density. 
This  leads  to  a  positive  value  of  yk.  Here,  we  use  the 
simplest  form, 

_  max(7"  )-min(T  ) 

yk~  y  = - 77U7 - ■  ,c\  .  (16) 

K  max(5i.)  —  minpt) 


to  illustrate  the  basic  methodology  of  this  analytical  ad¬ 
justment  procedure.  This  ratio  may  vary  with  depth.  A 
large  part  of  the  paper  by  Jackett  and  McDougall  (1995) 
was  devoted  to  developing  a  method  to  determine  yk. 
Interested  readers  are  referred  to  their  paper. 

Equations  (10),  (11),  (15),  and  (8)  represent  a  set  of 
2 K  algebraic  equations  for  2 K  unknowns  (A Tk,  ASk), 
k  =  1, 2, . . . ,  K.  Thus,  they  are  closure.  Among  them,  (8) 
is  nonlinear  and  (10),  (11),  and  (15)  are  linear. 

5.  Example 

The  same  example  as  described  in  section  2  is  used. 
Substitution  of  {Sk,  Tk.  zk}  values  from  Table  1  into  (8), 
(10),  (11),  and  (15),  and  the  Newton  iteration  method 
(Kelley  1987,  see  appendix  B)  is  used  to  solve  the  set  of 
2 K  algebraic  equations.  For  the  hydrographic  cast  listed 
in  Table  1,  only  three  iterations  are  needed  to  eliminate 
the  static  instability.  Tables  3  and  4  list  the  values  of 
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Table  4.  Change  of  ( Sk  +  (ppt)  at  each  iteration  using  Newton's  method.  It  is  noted  that  the  iteration  converges  at 

the  third  iteration. 


Depth  (m) 


7  =  0 


7  =  1 


7  =  2 


7  =  3 


7  =  4 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 
19 


0 

10 

20 

30 

50 

75 

100 

125 

150 

200 

250 

300 

400 

500 

600 

700 

800 

900 

1000 


34.424  300 
34.427  800 
34.288  000 
34.291  400 
34.299  100 
34.307  300 
34.328  000 
34.360  400 
34.369  600 
34.336  400 
34.341  500 
34.336  700 
34.285  200 
34.312  300 
34.402  200 
34.486  800 
34.490  400 
34.455  800 
34.475  500 


34.421  995 
34.420  749 
34.299  459 
34.298  031 
34.303  105 
34.309  152 
34.327  896 
34.358  223 
34.364  978 
34.335  234 
34.340  005 
34.333  394 
34.276  792 
34.320  875 
34.406  412 
34.488  540 
34.490  386 
34.454  421 
34.472  906 


34.421  985 
34.420  765 
34.299  526 
34.298  045 
34.303  102 
34.309  141 
34.327  883 
34.358  213 
34.364  978 
34.335  222 
34.339  993 
34.333  388 
34.276  820 
34.320  904 
34.406  410 
34.488  530 
34.490  374 
34.454  409 
34.472  897 


34.421  985 
34.420  765 
34.299  526 
34.298  045 
34.303  102 
34.309  141 
34.327  883 
34.358  213 
34.364  978 
34.335  222 
34.339  993 
34.333  388 
34.276  820 
34.320  904 
34.406  410 
34.488  530 
34.490  374 
34.454  409 
34.472  897 


34.421  985 
34.420  765 
34.299  526 
34.298  045 
34.303  102 
34.309  141 
34.327  883 
34.358  213 
34.364  978 
34.335  222 
34.339  993 
34.333  388 
34.276  820 
34.320  904 
34.406  410 
34.488  530 
34.490  374 
34.454  409 
34.472  897 


{Tk,  Sk }  at  the  each  iteration.  They  show  the  high  effi¬ 
ciency  of  this  method  for  elimination  of  static  instability 
in  the  hydrographic  cast.  Figure  2  shows  the  original  and 
adjusted  profiles 


{Sk,Tk,Ek},  k  =  l,2,...,K. 


(17) 


The  heat  and  salt  are  conserved  for  the  whole  water 
column  with  the  relative  root-mean  adjustment 

RRMA  =  0.0482.  (18) 

Comparing  (18)  to  (4),  we  may  find  that  this  analytical 
conserved  adjustment  scheme  has  a  smaller  RRMA 
(0.0482)  than  the  MA  method  (0.0712). 


Temp  difference:  -0.0000 


Salt  difference:  0.0000 


N2  difference:  -0.0000 


34.3  34.35  34.4  34.45  34.5 
Salinity 

Fig.  2.  As  in  Fig.  1,  but  using  the  analytical  conserved  method  proposed  in  this  paper. 
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Longitude 

Fig.  3.  Distribution  of  statically  unstable  casts  in  the  JPL-ECCO  10-day  data  centered  on  31  Dec  2008  (available  online  at  http://ecco.jpl. 
nasa.gov/external/).  The  data  were  produced  by  a  data  assimilation  system. 


Data  assimilation  is  required  in  operational  ocean 
data  access  and  retrieval  (Sun  1999).  It  is  to  blend  the 
modeled  variable  x„,  with  observational  data  yD  (e.g., 
Kalnay  2003;  Chu  et  al.  2004), 

+  W.[yo-//(Xffl)],  (19) 

where  x„  is  the  assimilated  variable,  H  is  an  operator 
that  provides  the  model’s  theoretical  estimate  of  what 
is  observed  at  the  observational  points,  and  W  is  the 
weight  matrix.  The  difference  among  various  data  as¬ 
similation  schemes  such  as  optimal  interpolation  (e.g., 
Lozano  et  al.  1996),  Kalman  filter  (e.g.,  Galanis  et  al. 
2006),  and  variation  methods  (e.g.,  Tang  and  Kleeman 
2004),  is  the  different  ways  to  determine  the  weight 
matrix  W.  The  data  assimilation  process  (19)  can  be 
considered  as  the  average  (in  a  generalized  sense)  of  x,„ 
and  yD.  In  ocean  ( T ,  5)  data  assimilation,  the  observa¬ 
tional  data  y0  may  contain  several  casts,  which  are 
statically  stable.  The  model  profile  x,„  is  also  statically 
stable  because  convective  adjustment  (Bryan  1969)  is 
usually  conducted  at  each  time  step. 

False  static  stability  may  be  generated  after  ( T ,  S )  data 
assimilation  [i.e.,  performing  (19)].  For  example,  10-day 
Jet  Propulsion  Laboratory  (JPL)  Estimating  the  Circu¬ 
lation  and  Climate  of  the  Ocean  (ECCO)  (7)  5)  fields 
centered  on  31  December  2008  (available  online  at 
http://ecco.jpl.nasa.gov/external/)  show  that  a  consider¬ 
able  portion  (35.32%)  of  profiles  are  statically  unstable 


(Fig.  3).  Here,  the  NODC’s  criterion  for  flagging  out 
statically  unstable  profiles, 

1—0.03  kgnT3  (0  >  zk  ^  —30  m) 

£min  =  |  -°-02  ks  m“3  (—30  m  >  Zk  S  -400  m) ,  (20) 

(o  kgnT'  (— 400m>zA.) 

is  used.  Because  such  a  false  static  instability  is  due  to 
the  blending  of  observational  data  with  the  model  data, 
it  is  not  a  real  instability.  Use  of  the  convective  adjust¬ 
ment  scheme  may  overcorrect  the  profiles. 

To  illustrate  this,  we  discuss  the  existing  convective 
adjustment  schemes  in  ocean  models.  The  various  con¬ 
vective  adjustment  schemes  are  based  on  the  same 
original  idea  (e.g.,  Bryan  1969):  whenever  a  water  col¬ 
umn  is  statically  unstable,  temperature  and  salinity  are 
vertically  adjusted  to  make  the  water  column  neutrally 
stable,  with  heat  and  salt  conserved  in  the  process.  The 
adjustment  takes  an  iterative  approach.  The  iteration 
continues  between  all  adjacent  levels  until  the  static 
instability  is  removed  in  the  whole  water  column.  Be¬ 
cause  the  adjustment  acts  on  only  neighboring  points, 
the  number  of  iterations  required  to  reach  the  final 
stable  state  is  infinite  for  a  given  unstable  profile  (Smith 
1989).  In  practice,  however,  the  number  of  iteration  is 
always  finite,  and  this  leads  to  some  residual  instability 
(Killworth  1989). 
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Temp  difference:  0.0000  Salt  difference:  0.0000  N2  difference:  0.0773 


Fig.  4.  As  in  Fig.  1,  but  using  the  complete  convective  adjustment  method  (Yin  and  Sarachik  1994). 


Several  algorithms  were  developed  to  remove  these 
residual  static  instabilities  such  as  the  implicit  vertical 
diffusion  scheme  (Cox  1984;  Killworth  1989)  and  the 
complete  adjustment  scheme  (Yin  and  Sarachik  1994). 
The  former  tests  the  static  stability  between  the  verti¬ 
cally  adjacent  levels  and,  if  unstable,  the  vertical  diffu- 
sivity  is  set  to  a  large  value  (convective  diffusivity)  to 
smooth  out  the  instability.  The  latter  determines  the 
upper  and  lower  boundaries  of  each  adjusted  region 
while  keeping  the  instantaneous  adjustment  within  each 
unstable  region.  Yin  and  Sarachik  (1994)  showed  that 
the  complete  convective  adjustment  scheme  is  more 
efficient  than  the  implicit  vertical  diffusion  scheme  and 
guaranteed  a  complete  removal  of  static  instability  of  a 
water  column  at  each  time  step.  For  the  same  example  as 
described  in  section  2,  the  complete  convective  adjust¬ 
ment  scheme  removes  the  static  instabilities  (Fig.  4)  with 
the  relative  root-mean  adjustment 

RRMA  =  0.2192.  (21) 

This  value  is  4.5  times  larger  than  that  of  (0.0482)  using 
the  analytical  adjustment  method. 

6.  Conclusions 

A  new  analytical  conserved  adjustment  scheme  is 
developed  to  eliminate  the  static  instability  of  raw  and 
averaged  observational  hydrographic  data.  This  method 


adjusts  the  temperature  and  salinity  profiles  (A  Tk,  A  .S'/,, 
k  =  1,  2,  . . .  ,  K}  simultaneously  and  efficiently  on  the 
basis  of  three  types  of  constraints:  1)  heat  and  salt  con¬ 
servation,  2)  predetermined  ( ATkIASk )  ratios  (or  called 
adjustment  ratios)  for  all  levels,  and  3)  the  removal  of 
static  instability  by  adjusting  the  static  stability  with 
a  combined  conservation  and  nonuniform  increment 
treatment.  With  these  constraints,  a  set  of  2 K  combined 
linear/nonlinear  algebraic  equations  are  established  for 
(A Tk,  A. S’/.}.  Among  them,  (K  +  i )  algebraic  equations 
are  linear,  and  ( K  —  1)  equations  are  nonlinear.  Newton’s 
method  is  used  to  solve  this  set  of  equations.  The  pro¬ 
posed  scheme  has  very  small  relative  root-mean-square 
adjustment  compared  to  the  existing  methods.  More¬ 
over,  it  has  three  features:  1)  conservation  of  heat  and 
salt,  2)  removal  of  static  instabilities  with  small  ( T ,  S) 
adjustments,  and  3)  analytical  form.  With  these  features, 
it  can  be  widely  used  in  ocean  ( T ,  S)  data  analysis.  Be¬ 
sides,  ocean  data  assimilation  may  cause  false  static 
instabilities.  Because  this  instability  is  not  real,  com¬ 
monly  used  convective  adjustment  schemes  may  over¬ 
adjust  the  profiles.  Therefore,  the  proposed  analytical 
conserved  scheme  can  be  used  in  ocean  (T,  S)  data 
assimilation. 
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APPENDIX  A 

Validity  of  the  Conservation  Constraints 

In  ocean  modeling,  all  the  convective  adjustment 
schemes  for  stabilizing  ( T ,  5)  profiles  require  heat  and 
salt  conservation  (e.g.,  Yin  and  Sarachik  1994).  In  ocean 
data  analysis,  such  conservation  constraints  are  also 
valid.  After  quality  control  procedures,  it  is  reasonable 
to  assume  that  ocean  observational  data  i (i  contain  ran¬ 
dom  error  ijj', 


individual  observations  <xe.  Thus,  the  conservation  con¬ 
straints  (10)  and  (11)  guarantee  that 

0Aadj)  =  0‘)  +  m,  (as) 

the  same  smaller  error  variance  of  the  vertically  in¬ 
tegrated  observed  values. 

APPENDIX  B 


(Al) 


Newton  Method 


with  population  mean  (i //')  =  0  and  standard  deviation 
cre.  Here,  ip  denotes  ( T ,  .S’),  and  if/  is  the  true  value  at  the 
same  location  and  time  as  the  observation  taken  place. 
The  population  mean  of  (Al)  gives 

(*)  =  W).  (A2) 

An  observational  profile  (ipk,  k  =  1,  2,  . . .  ,  K)  can  be 
taken  as  a  sample.  Vertical  integration  of  the  observa¬ 
tional  profile  is  represented  by  weighted  average  [see 
(12)], 

K  K  K 

(<A)  =  0)  =  (t/0  =  Y^akijj'k. 

k=  1  k= 1  k=  1 

(A3) 

The  random  errors  at  different  depth  i/4  are  considered 
independent.  The  central  limit  theorem  states  that  the 
linear  combination 


Let  the  temperature  and  salinity  adjustment  be  rep¬ 
resented  by  a  2K-dimensional  vector  P, 


’  Pi 

-A7V 

Pi 

ASi 

Pi 

^ t2 

Pa 

= 

\s2 

Pm- i 

^tk 

.  Pm  . 

1 

> 

M  =  2K. 


(Bl) 


The  algebraic  Eqs.  (10),  (11),  (15),  and  (8)  [note  that  we 
put  (8)  at  the  last]  can  be  represented  by 

F(P)  =  0,  (B2) 


K 


Y'  = 


where  F  has  the  dimension  of  2 K.  The  classical  Newton 
method  (Kelley  1987)  for  approximating  a  desired  so¬ 
lution  P  to  (B2)  is  formally  defined  by  the  iteration 


has  a  normal  distribution  with  zero  mean  and  variance,  pO+O  =  p(;)  _  j  1(p0'))p(p(;)),  j  =  0,  1,  2,  ...  ,  (B3) 


K  K 

<r Y'  =  ^aWe  =  (A5) 

k= 1  k= 1 

From  (14),  we  have 

K  K 

YjoI  <  ^hak  =  1.  (A6) 

k= 1  k= 1 

Substitution  of  (A6)  into  (A5)  leads  to 

<t  K,<  (j  e,  (A7) 

which  indicates  that  the  error  variance  of  the  vertically 
integrated  observed  values  i/f  is  smaller  than  that  of  the 


where  P(,)  is  the  ;th  approximation  to  the  solution  of 
(B2),  J(P*A)  is  the  Jacobian  matrix  of  F(P)  evaluated  at 
P(y\  Inversion  of  the  Jacobian  matrix  is  not  performed 
in  practice;  rather  (B3)  is  implemented  via  solution  of 
the  following  system  of  linear  equations  at  the  each 
iteration: 

j(p(/))  .  d(h  =  b0),  b(y)  =  — F(P(y)),  (B4) 

followed  by  the  update 

p(/+i)  =  pO)  +  dW,  (B5) 

where  d(;)  is  called  the  Newton  direction.  The  iteration 
stops  at  step  /  when 
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max|<i^|  <  10  6  (K  for  temperature 

and  ppt  for  salinity).  (B6) 

When  the  set  of  algebraic  equations  take  the  order  of 
(10),  (11),  (15),  and  (8),  the  Jacobian  matrix  J(P^)  with 
dimension  of  M  X  M  is  represented  by 


■flu 

a\2 

a\M 

J(P0))  = 

U21 

a22 

a2M 

_aMl 

aM2  ' 

■  aMM . 

J(PW) 


* 

0 

* 

0 

* 

0 

* 

0  ■ 

0 

* 

0 

* 

0 

* 

.  0 

* 

* 

* 

0 

0 

0 

0 

.  0 

0 

* 

* 

* 

* 

0 

0 

.  0 

0 

0 

0 

* 

* 

0 

0 

.  0 

0 

0 

0 

* 

* 

* 

* 

.  0 

0 

0 

0 

0 

0 

*  *  0 

0 

0 

0 

0 

0 

*  *  * 

* 

(B8) 


where  the  M  X  M  elements  are  given  in  (Cl)  of  ap¬ 
pendix  C.  The  Jacobian  matrix  (B7)  has  the  following 
format  with  many  zero  elements: 


where  nonzero  elements  are  indicated  by  asterisks.  The 
vector  b  in  the  right-hand  side  of  (B4)  has  the  following 
components: 


bx  =0,  b2  =  0,  b3  =0,  b4=  E**  -  p(5j  +  A S[»,  T]  +  A +  P(S2  +  A S">,  T2  +  A 7™,  z2), 
b5  =  0,  b6  =  E '**  -  P(S2  +  A T2  +  A  lf,z2)  +  P(S3  +  A  $>,  T3  +  A  Tf,z2),  bM_3  =  0, 

bM  =  E**,-p{SK_,  +  *!$-v  Tk  +  A  l*£-vZK-d  +  P(SK  +  AS«,  T  K  +  A 


r(J) 


c(j) 


rU) 


(B9) 


It  is  noted  that  the  well-posed  linear  algebraic  Eq.  (A4) 
is  easily  solved  with  the  initial  guess, 

P(°)  =  0,  (B10) 

that  is, 

Ar[0)=0,  A5j°}  =  0,  k=l,2,  ...  ,K.  (Bll) 

With  the  initial  guess  (Bll),  the  Newton  direction  d<()) 
is  obtained  from  solving  the  linear  algebraic  Eq.  (B4). 
The  vector  d(0>  is  added  to  the  initial  guess  P(0\  which 
leads  to 


the  adjustment  stops.  Otherwise,  the  iteration  continues, 
that  is,  the  linear  algebraic  Eq.  (B4)  is  solved  after  using 
pB)  from  (B12).  Addition  of  the  solution  d*1-1  to  P(1) 
leads  to  P(2).  If  there  is  no  static  instability,  the  adjust¬ 
ment  stops.  Otherwise,  the  iteration  continues  until  the 
static  instability  is  eliminated. 


APPENDIX  C 

Elements  of  Jacobian  Matrix  (B7) 

The  elements  of  M  X  M  Jacobian  matrix  (B7)  are 
given  by 


P(t)  =  P(°)  +  d(«)  (B12) 


With  P(1\  the  cast  is  adjusted  to  its  first  iterated 
values. 


S^  =  Sk  +  AS^\  if  =  Tk  +  ATjj?.  (B13) 

Substitution  of  (B12)  into  (1)  gives  static  stability  after 
the  first  iteration  Ek  \  If 

*=1,2. . K,  (B14) 


A  z,  A  z,  +  Az7 

Cl  —  - —  Cl  —  - - - — 

11  2  ’  13  2  ’  ■ "  ’ 

_  ^N- 2  +  ^N-l  „  _  A ZN_t 

2  ’  UhM-l  2  ’ 

fl12  =  fl14  =  fl16  -  '  '  '  =  a\M  =  0’ 


A  z,  A  z,  +  A  z2 

n  =  - L  n  =  - i - £ 

22  2  ’  24  2 


J  ■  •  •  ) 


A^at-2  +  A  zN x 


” 2,M-2  2  ’  2 

fl21  =  a23  =  fl25  =  ''  '  =  a2,M-l  = 


AZiV l 


fl31  fl32  fl33  fl34  fl35  "  fl3  M 
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dp(S1+AS<^\T1  +  T(/),z2) 
dT 

- P(S !  +  A S[j\  T1  +  A T[i),z2)a(S1  +  A s[D, 

T,+AT[»,z2), 

dP(S j  +  AS<y),  7\  +  A7f\z2) 
dS 

P(S1  +  A S<y),  Tj  +  A T[i),z2MS1  +  AS[n, 

Tl+AT^,z2), 

dp(S2  +  AS^,T2  +  AT^,z2) 
dT 

p(S2  +  A  S(J\  T2  +  A  lf,z2MS2  +  AA<y), 

T2  +  A  T^,  z2), 

dp(S2  +  AS^,T2+AT^\z2) 

9S 

-p(S2  +  AS2  \  T2  +  AT2\z2)f3(S2  +  A sf, 
T2  +  ATf,z2), 

a46=  "  '  =  a4M  = 

1,  fl54  —  72’ 

fl52  =  °55  =  a56  =  fl57  =  '  '  '  =  = 

dp(S2  +  AS^,T2  +  ATf,z3) 
dT 

-p(S2  +  AS2  \  T2  +  A if,  z3)a{S2  +  A S(J\ 
T2  +  AT(J\z3), 


aM-l,M-3  ~  aM-\,M-2  ~  ^N-V 

flM-l,l  —  aM- 1,2  —  ■  ■  ■  —  aM_i  M_ 4 
~  aM—l,M—l  ~  aM~l,M  ~ 

_  dp(SK _j  +  A +  A T{J)_1,zk) 

aM,M- 3  gj 

=  -p(Vi  +  AS{lv  Tk_3  +  AT{lvzK)a(SK_, 

+  AS(lvTK_l  +  AT^_vzK), 

_  dP(SK_ +  A.S«1;  7-^,  +  A T^_vzk) 

UM,M- 2  35 

=  p(V,  +  AS«  lf  Tk1  +  AT«  1>Zjf)j8(SJC_1 

+  A5y)_1,rx_1  +  Arw_1,zx), 

ap^  +  As^.r^  +  Ar^.Zjc) 

aM.M- 1  g'p 

=  P(S*  +  AS«  7*  +  ATf,zK)a(SK  +  A  S$, 
Tk  +  M$,zk), 

dp(SK  +  AS<i\TK  +  AT(j?,zK) 

aMM  gS 

=  -P(SK  +  A S$,  Tk  +  AT$,zk)KSk  +  AfiW, 
tk  +  Atk^k),  and 


fl2.M  '  '  '  aM-4,M 


(Cl) 
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